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Abstract 

RNA interference (RNAi) is a eukaryotic molecular system that serves two primary functions: 1) gene regulation and 2) protection 
against selfish elements such as viruses and transposable DNA. Although the biochemistry of RNAi has been detailed in model 
organisms, very little is known about the broad-scale patterns and forces that have shaped RNAi evolution. Here, we provide a 
comprehensive evolutionary analysis of the Dicer protein family, which carries out the initial RNA recognition and processing 
steps in the RNAi pathway. We show that Dicer genes duplicated and diversified independently in early animal and plant 
evolution, coincident with the origins of multicellularity. We identify a strong signature of long-term protein-coding adaptation 
that has continually reshaped the RNA-binding pocket of the plant Dicer responsible for antiviral immunity, suggesting an 
evolutionary arms race with viral factors. We also identify key changes in Dicer domain architecture and sequence leading to 
specialization in either gene-regulatory or protective functions in animal and plant paralogs. As a whole, these results reveal a 
dynamic picture in which the evolution of Dicer function has driven elaboration of parallel RNAi functional pathways in animals 
and plants. 

Key words: RNAi, RNA interference, dicer, positive selection, antiviral dicer, PAZ domain. 



Introduction 

RNA interference (RNAi) is a widespread molecular process 
that plays two primary cellular roles: 1) regulating the activ- 
ities of endogenous genes and 2) protecting the cell from 
selfish genetic material such as viruses and transposable elem- 
ents. RNAi's gene-regulatory function is largely facilitated by 
degradation of cytoplasmic mRNAs by the RNA-induced 
silencing complex (RISC), a collection of proteins and tem- 
plate RNA that targets and degrades complementary RNA 
sequences (Zamore et al. 2000; Preall and Sontheimer 2005; 
Liu et al. 2006; Wang, Noland, et al. 2009). RNAi is additionally 
thought to affect gene regulation by altering chromatin struc- 
ture (Hall et al. 2003; Volpe et al. 2003; Fukagawa et al. 2004). 
The protective role of RNAi appears to occur primarily 
through direct identification, processing and degradation of 
unwanted RNA molecules by the RISC (Galiana-Arnoux et al. 
2006). 

The RNAi process occurs through an interaction between 
RNA molecules and multiple proteins that collectively 
identify, process, and regulate RNA. Dicer is the primary 
RNA recognition and processing protein. Dicer anchors a 
double-stranded RNA (dsRNA) molecule and cuts it to pro- 
duce short dsRN As (Ketting et al. 2001 ). These short RN As are 
converted to single-stranded form and bound to an 
Argonaute (AGO) protein through a process coordinated 
by Dicer and other RNA-binding proteins (Hammond et al. 



2001). Once bound to AGO, the template RNA binds com- 
plementary RNA sequences, at which point they are cleaved 
or otherwise regulated by AGO and other components of 
the RISC (Wang, Juranek et al. 2009). 

Because of its primary role in dsRNA recognition and pro- 
cessing, Dicer has received considerable research attention 
(Merritt et al. 2010; Rossbach 2010; Pecot et al. 2011; 
Davalos and Esteller 2012; Grimm 2012). The typical Dicer 
found in animals and plants is a large protein consisting of 
a series of functional domains that interact with RNA in 
different ways (fig. 1). The PAZ domain anchors the 3'-end 
of the dsRNA, after which it is cleaved by twin RNase3 
domains (Bernstein et al. 2001). This "catalytic core" is suffi- 
cient for Dicer function in some organisms (Macrae et al. 
2006). Some Dicer proteins additionally anchor the 5' end 
of the dsRNA molecule using an N-terminal extension 
of the PAZ domain. This S f anchoring seems particularly 
important for micro-RNA (miRNA) recognition and process- 
ing, supporting Dicer's main gene-regulatory role (Park 
et al. 2011). 

Dicer's DEAD/Helicase domain is thought to facilitate 
the movement of the protein along long dsRNA molecules 
(Welker et al. 2011), suggesting it may be important for sup- 
porting Dicer's protective role by allowing it to find the end 
of viral and transposable-element RNAs, which can then be 
assayed by the PAZ domain. Further support for this model 
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Fig. 1. Dicer protein consists of multiple RNA-interacting functional domains. We plot the domain architecture of human Dicer along the protein 
sequence; other Dicer proteins from animals and plants have similar domain architectures (figs. 2 and 4). Functional domains were identified via 
sequence searches of the PFam database (Punta et al. 2011), the SMART database (Letunic et al. 2012) and the NCBI conserved-domain database 
(March ler-Bauer et al. 2011). The 5 r RNA-binding pocket extension of the canonical PAZ domain found in some animal Dicers is also shown (Park et al. 
2011). 



comes from examination of Drosophila melanogaster Dicerl, 
which seems to have specialized in miRNA processing by 
losing its functional DEAD/Helicase domain (Welker et al. 
2011). Other Dicer functional domains appear to coordinate 
the "hand-off" of processed RNAs to AGO, either through 
direct Dicer-RNA interaction or through interactions with 
other partner proteins (Maniataki and Mourelatos 2005; 
Koscianska et al. 2011). 

Although the biochemical functions of Dicer have been 
detailed in model organisms, the evolution of the Dicer 
"superfamily" remains poorly characterized. Dicer is absent 
from bacteria and archaea but is found throughout eukary- 
otes, suggesting an early eukaryote origin (Cerutti and 
Casas-Mollano 2006; Shabalina and Koonin 2008). Current 
evidence suggests that the Dicer family diversified independ- 
ently in animals, plants, and fungi (Cerutti and Casas-Mollano 
2006) and was lost from many parasitic protozoa (Ullu et al. 
2004; Baum et al. 2009) as well as model fungi lacking RNAi 
(Drinnenberg et al. 2009). However, the support in favor of 
this model is relatively weak, and alternative hypotheses have 
not been thoroughly evaluated. 

Vertebrates and nematodes have only one Dicer gene, 
whereas insects have two (Hammond 2005), suggesting an 
insect-specific duplication followed by functional divergence 
into miRNA-based gene regulation and antiviral immunity 
(de Jong et al. 2009). This hypothesis is supported by evidence 
for strong positive selection affecting fly Dicer2 — which per- 
forms an antiviral function (Obbard et al. 2006; Heger and 
Ponting 2007; Kolaczkowski et al. 201 1) — and a parallel loss of 
DEAD/Helicase function in Dicerl, which appears to focus 
this protein's function on miRNA processing (Welker et al. 
2011). All of this is consistent with a model of gene duplica- 
tion followed by functional divergence in insects or arthro- 
pods. However, phylogenetic analysis — the real test of 
macro-evolutionary hypotheses (Huelsenbeck and Rannala 
1997) — has so far failed to strongly support the insect-specific 
duplication hypothesis (de Jong et al. 2009). 

Most model plant genomes encode four Dicer genes (DCLs 
1-4), which — similar to the case in animals — appear to have 
diverged to function in miRNA-based gene regulation vs. anti- 
viral immunity (Blevins et al. 2006; Bouche et al. 2006). 
However, there may be some functional overlap among 
plant Dicer paralogs, particularly in the case of antiviral 
Dicers, where one Dicer may compensate for loss of a para- 
log's function (Gasciolli et al. 2005). How plant Dicers func- 
tionally diverged is completely unknown, so it is impossible to 
evaluate whether there is any similarity with what we observe 
in animals. 



Here, we examine the broad patterns of Dicer evolution 
using a combination of phylogenetic, structural-modeling and 
sequence-analysis approaches. We show that: 1) Dicer inde- 
pendently diversified in animal and plant lineages, coincident 
with the origins of multicellularity and requirements for com- 
plex gene regulation; 2) animal Dicer did not duplicate in 
insects but much earlier in metazoan evolution, with antiviral 
Dicer2 being subsequently lost from lineages developing al- 
ternative antiviral strategies; 3) the main plant antiviral Dicer 
(DCL-4) has been a repeated target of intense positive selec- 
tion for changes in RNA recognition and/or binding, suggest- 
ing a long-term evolutionary arms race between this protein 
and viral molecules; and 4) although the biochemical capacity 
to recognize miRNAs appears ancestral, efficient miRNA rec- 
ognition like that employed by humans arose later and pos- 
sibly independently in animals and plants. These results 
provide a thorough picture of the forces and patterns shaping 
Dicer evolution and suggest that many common assumptions 
about the evolution of RNAi may warrant more careful 
investigation. 

Results 

Evolution of Eukaryote Dicers 

The availability of complete genome sequences from a variety 
of eukaryotes provides an opportunity to systematically 
examine the origin and evolution of important biomolecular 
systems by studying the evolution of the genes encoding their 
component parts. To deepen our understanding of how the 
RNAi system evolved, we examined the evolution of the Dicer 
gene family, which carries out the initial RNA recognition and 
processing steps that initiate RNAi (Hamilton and Baulcombe 
1999; Bernstein et al. 2001). 

We used homology-based gene identification based on 
confirmed functional Dicers to find novel Dicer genes from 
a large number of fully sequenced genomes (supplementary 
Excel spreadsheet, Supplementary Material online). As ex- 
pected, we found full-length Dicer genes throughout animals, 
plants, fungi, and many protozoan lineages, although Dicers 
were absent from yeast and some parasitic protozoa, consist- 
ent with previous results (Cerutti and Casas-Mollano 2006). 
The functional domain building-blocks of Dicer proteins are 
encoded in bacteria and archaea, but we found no full-length 
or multi-domain Dicer proteins in these lineages, suggest- 
ing that Dicer — and indeed the entire canonical RNAi 
machinery — is a eukaryote innovation that arose very early 
in this lineage's evolution. 

We inferred the evolutionary history of eukaryote Dicers 
using maximum-likelihood and Bayesian techniques, finding 
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Fig. 2. Phylogenetic analysis supports independent expansions of an 
ancient eukaryote Dicer protein in animals and plants. We plot the 
support for monophyletic expansions of Dicer paralogs in animals 
and plants. Support is given as SH-like aLRT scores/maximum-likelihood 
bootstrap proportions/Bayesian posterior probabilities. See supplemen- 
tary figures S1-S5, Supplementary Material online, for full trees and 
additional analyses. 

strong support for independent expansions in animals, plants, 
and fungi (fig. 2 and supplementary fig. S1, Supplementary 
Material online). These independent expansions were consist- 
ently supported by multiple statistical-support measures 
and alignment strategies, suggesting our results are robust 
to alignment ambiguity and do not depend on any particular 
approach for inferring clade support. Support for key 
nodes increased when we used an elision strategy that com- 
bines information from multiple alignments (Wheeler et al. 
1995), suggesting that the most reliable alignment positions 
consistently support a phylogeny in which Dicer genes inde- 
pendently expanded in the major multicellular eukaryote 
lineages. 

Support for independent Dicer expansions increased when 
we removed the fast-evolving taxa most likely to contribute 
to phylogenetic error, suggesting our result is robust to taxon 
sampling and not strongly affected by long-branch attraction 
(supplementary fig. S2, Supplementary Material online). 
Our phylogeny was also robust to evolutionary model 
uncertainty. Alternative evolutionary models identified 
using model-selection procedures all produced the same 
maximum-likelihood phylogeny (supplementary table S1, 
Supplementary Material online), and Bayesian analyses that 
integrate over multiple evolutionary models recovered the 
same independent expansions of animal, plant, and fungi 
Dicers with high posterior probability (supplementary fig. 
S1, Supplementary Material online). 

We also reconstructed the Dicer gene family tree using 
advanced heterogeneous evolutionary models that have 
been shown to prevent errors and improve phylogenetic ac- 
curacy under realistic conditions. First, we used a Bayesian 
technique that incorporates heterogeneity in amino acid fre- 
quencies across the protein sequence (Lartillot et al. 2009). 
This approach recovered the same independent Dicer expan- 
sions in animals, plants, and fungi, although statistical support 



was somewhat reduced, as expected with a more complex 
model (supplementary fig. S3, Supplementary Material 
online). Animal and plant expansions were also present in 
the maximum-likelihood phylogeny reconstructed using a 
model that incorporates across-site branch-length heterogen- 
eity (Kolaczkowski and Thornton 2008), although fungi Dicers 
were not monophyletic in this tree (supplementary fig. S4, 
Supplementary Material online). 

We were concerned that the high variability of the overall 
Dicer domain architecture across eukaryotes (supplementary 
Excel spreadsheet, Supplementary Material online) could lead 
to phylogenetic artifacts, so we additionally reconstructed 
the Dicer phylogeny using only the aligned "catalytic core," 
consisting of the PAZ and twin RNase3 domains conserved 
across all Dicers. This phylogeny also supported independent 
Dicer expansions in animals, plants, and fungi with high 
statistical confidence (supplementary fig. S5, Supplementary 
Material online). 

Together these results strongly support an evolutionary 
history in which the Dicer gene family originated early in 
eukaryote evolution and independently expanded in ani- 
mals, plants and fungi, coincident with the evolution of multi- 
cellularity — and requirement for more complex gene 
regulation — in these major eukaryote lineages. This result 
is robust to all the major factors that have been shown to 
affect phylogenetic accuracy, including the reconstruction 
strategy (maximum-likelihood vs. Bayesian), the method 
for assessing statistical confidence, alignment ambiguity, 
taxon sampling and evolutionary model ambiguity and 
inadequacy. 

Our overall Dicer phylogeny provides some clues about the 
patterns of expansions in animals and plants (supplementary 
figs. S1-S5, Supplementary Material online). Specifically, we 
observed a main clade of animal Dicerl homologs spanning 
arthropods, nematodes, and vertebrates as well as a smaller 
clade of Dicer2 homologs from arthropods, which fell toward 
the unresolved base of the animal tree. In plants, the main 
DCL types (DCLs 1-4) appear to have originated early, as 
monocots and dicots are represented in all four DCL clades, 
and DCLs 1 and 3 contain clear moss and Selaginella 
orthologs. 

However, many of the fine-scale branching patterns 
remained unresolved in the full Dicer tree, limiting our ability 
to draw strong conclusions about the patterns of Dicer evo- 
lution within animals and plants. We hypothesized that align- 
ing highly divergent Dicer homologs from across eukaryotes 
reduced the alignable data set size, which could contribute to 
lack of resolution in the full tree. Indeed, after processing 
alignments to remove ambiguously aligned regions, the com- 
bined sequence alignment was much shorter than individual 
alignments of only animal or plant sequences, and clade 
support was noticeably lower (supplementary table S2, 
Supplementary Material online). Because our full Dicer tree 
strongly supported independent expansions in animals and 
plants, we reconstructed separate animal and plant phyloge- 
nies to make better use of the available sequence data 
and generate a clearer picture of precisely when and how 
these independent expansions occurred. 
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Evolution of Animal Dicers 

We reconstructed the evolutionary history of animal Dicers 
using maximum-likelihood and Bayesian methods (fig. 3 
and supplementary fig. S6, Supplementary Material online). 
We found strong support for an ancient duplication, giving 
rise to Dicerl and Dicer2 genes very early in metazoan evo- 
lution. This is in stark contrast to the widely held belief that 
the Dicerl /2 duplication occurred later in insects (de Jong 
et al. 2009). 

As with our evolutionary analysis of eukaryote Dicers as a 
whole, we evaluated the robustness of our animal Dicer tree 
to all the major factors known to affect phylogenetic accur- 
acy. We found that our results were robust to alignment 
ambiguity and use of different methods for assessing statis- 
tical confidence (supplementary fig. S6, Supplementary 
Material online). Support for an ancient Dicer duplication 
in animals remained after removing fast-evolving taxa (sup- 
plementary fig. S7, Supplementary Material online) and 
was robust to model uncertainty (supplementary table S1 
and fig. S6, Supplementary Material online). Use of advanced 
heterogeneous evolutionary models less susceptible to phylo- 
genetic errors also recovered this ancient Dicerl /2 duplication 
event (supplementary figs. S8 and S9, Supplementary Material 
online), as did reconstructing the animal Dicer phylogeny 
using only the aligned "catalytic core" of the protein (supple- 
mentary fig. S10, Supplementary Material online). 

Although it is impossible to completely rule out phylogen- 
etic error in any analysis, we were unable to detect a strong 
affect on our results from any of the major factors known to 
cause such errors, suggesting our results are indeed reliable. 
Furthermore, we found that non-insect animals — the prawn 
Litopenaeus vannamei and the planarians Clonorchis sinensis 
and Schmidtea medeterranea — had two Dicer genes, one of 
which grouped with insect Dicerl, and the other of which 
strongly grouped with insect Dicer2, indicating that the 
Dicerl /Dicer2 duplication occurred prior to insect divergence 
(fig. 3). Statistical topology tests also strongly rejected an 
arthropod- or insect-specific gene duplication as a possible 
explanation for our data, using full-length Dicer alignments 
(P < 3.13 x 10~ 3 ), alignments of only the DEAD/Helicase do- 
mains (P < 2.65 x 10 -2 ) or alignments of the conserved Dicer 
"catalytic core" (P<1.01 x 10 -3 , supplementary table S3, 
Supplementary Material online). 

Synthesizing results from multiple phylogenetic analyses of 
various alignments, we identified Dicer2 orthologs from the 
basal metazoa, Trichoplax and Nematostella, suggesting that 
the Dicerl /2 duplication may have occurred very early in 
animal evolution. However, the relatively sparse taxon sam- 
pling and long-branch lengths at the base of the animal phyl- 
ogeny left much of the basal branching pattern unresolved. 

To attempt to better resolve this early branching pattern, 
we aligned additional animal Dicer RNase3 sequence data 
from a previous study (de Jong et al. 2009) to homologous 
regions from our sequence database and reconstructed the 
animal Dicer phylogeny from this expanded dataset. 
Although the same early Dicerl /2 duplication was recovered, 
the additional data did not add much taxonomic diversity, 



and the basal branching pattern remained weakly resolved 
(supplementary fig. S11, Supplementary Material online). 

In addition to retrieving data from additional taxa, phylo- 
genetic resolution can often be improved by including add- 
itional sequence data; unfortunately, this approach is not 
applicable to gene-tree reconstruction, as the length of the 
gene's sequence cannot be extended. However, it has long 
been observed that physically and functionally interacting 
proteins often have highly correlated patterns of gene dupli- 
cation (Fryxell 1996; Pazos et al. 1997; Goh et al. 2000; Wu 
et al. 2011). Based on this observation, if the phylogenies of 
genes that interact in the same biochemical pathway all ex- 
hibit the same duplication pattern, this not only suggests that 
the entire pathway duplicated at the same time but also 
reinforces our confidence that the individual gene trees are 
robust, as phylogenetic artifacts would have to strongly affect 
all the gene trees in the same way to recover a consistent — 
but erroneous — phylogeny. 

In the case of RNAi, Dicer's function is tightly linked to 
AGO proteins, with which Dicer directly interacts (van Rij 
et al. 2006; Marques et al. 2010). Furthermore, studies in 
model organisms have shown that Dicerl interacts with 
AGOl, whereas Dicer2 primarily interacts with AG02 (Liu 
et al. 2003; Okamura et al. 2004). When we reconstructed 
the AGO gene family tree by maximum likelihood, we found 
strong support for an AG01/AG02 gene duplication occur- 
ring at the same time in early animal evolution that we 
observed in our Dicer phylogeny (supplementary fig. S12, 
Supplementary Material online). We also observed a striking 
correlation between the taxonomic distributions of 
Dicer-AGO interacting partners, with Dicer2/AG02 found 
in arthropods, prawn, planaria, and Trichoplax. Our AGO 
tree also recovered a strongly supported Trichoplax AG02 
along with a strongly supported Nematostella AGOl, further 
suggesting that this duplication event may have occurred very 
early in animal evolution. This result reinforces our confidence 
in an ancient Dicerl /2 duplication and suggests that — in add- 
ition to Dicer — many of the core proteins involved in RNAi 
may have duplicated around the same time early in animal 
evolution. 

Under this early duplication scenario, both Dicer2 and 
AG02 appear to have been lost from deuterostomes and 
nematodes. These genes may have also been lost from mol- 
lusks/annelids, although more extensive taxon sampling from 
these lineages will be required to definitively address this issue. 
We also observed losses of Dicerl and AGOl from Trichoplax, 
suggesting that this simple organism has lost the RNAi ma- 
chinery required for miRNA-based gene regulation. These 
findings are consistent with recent studies suggesting 
Trichoplax lacks miRNA genes (Grimson et al. 2008; Hertel 
et al. 2009). 

It appears that the Dicerl /2 and AGOI/2 duplication 
events both occurred after sponges diverged from the main 
animal branch but before the cnidarian split. However, the 
precise timings of the Dicerl /2 and AGOI/2 duplications are 
difficult to confidently ascertain, given the long branch 
lengths and lack of resolution at the base of the animal phyl- 
ogeny in both the Dicer and AGO trees. 
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Fig. 3. Dicer duplicated early in animal evolution. We inferred the metazoan Dicer family phytogeny using maximum likelihood and Bayesian methods. 
Support is shown for key nodes as SH-like aLRT scores/bootstrap proportions/Bayesian posterior probabilities (supplementary figs. S6-S11, 
Supplementary Material online, for additional support calculations and analyses). We also show the inferred domain architecture for each sequence. 
We used a branch-sites model to identify protein-coding adaptation from aligned codon sequences (Zhang et al. 2005). Colored circles indicate 
branches showing significant support for adaptation in each functional domain (P < 0.05 after correcting for multiple tests). Species names are colored 
by taxonomic group. 
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The domain architecture of large, multi-domain proteins 
can change over evolutionary timescales, leading to major 
changes in the protein's functional repertoire. To examine 
the evolution of Dicer's domain architecture in animals, 
we identified known functional domains in each sequence 
using the PFam database (Punta et al. 2011) and confirmed 
each inferred domain architecture using the SMART 
(Letunic et al. 2012) and CDD (March ler-Bauer et al. 
2011) databases (supplementary Excel spreadsheet, 
Supplementary Material online). As figure 3 shows, we 
observed a complete loss of the DEAD functional domain 
in arthropod, mollusk, annelid, and planarian Dicerl. 
Additionally, we observed a corruption of the Walker-B 
"DECH" motif in the majority of these sequences, which 
is likely to disrupt the ATP-binding properties of this 
domain and lead to a complete loss of DEAD/Helicase 
function (Sengoku et al. 2006). Nematode Dicerl retained 
a recognizable DEAD/Helicase domain that appears to be 
functional. Since the loss of a functional DEAD/Helicase 
domain in fly Dicerl has been linked to specialization 
in miRNA-based gene regulation (Welker et al. 2011), our 
results suggest that this functional specialization may have 
occurred multiple times in various protostome lineages, 
including early in arthropod evolution. 

Previous studies have found that antiviral Dicer2 is under 
intense positive selection in Drosophila melanogaster and 
across the Drosophila phylogeny (Obbard et al. 2006; Heger 
and Ponting 2007; Kolaczkowski et al. 2011). We confirmed 
that Dicer2 DEAD/Helicase and PAZ domains have experi- 
enced positive selection in flies using branch-sites analyses to 
identify adaptive protein-coding changes (Zhang et al. 2005). 
However, other than a spattering of adaptive substitutions 
in various lineages, we found little evidence for recurrent 
adaptation across the animal Dicer phylogeny (fig. 3). 
We note that adaptively driven loss of domain function is 
unlikely to be detectable using branch-sites models based 
on the nonsynonymous/synonymous substitution rate ratio 
(dN/dS), as any signal will be rapidly lost due to the fast 
evolution of nonfunctional protein sequences. For this 
reason, it is not possible to determine whether the loss of 
DEAD/Helicase function observed in various protostomes 
was adaptive. 

Evolution of Plant Dicers 

We used homology-based gene prediction to identify 
Dicer genes from fully sequenced plant genomes and recon- 
structed their phylogeny using maximum-likelihood and 
Bayesian methods. We found that the Dicer gene underwent 
a rapid four-way duplication early in plant evolution, giving 
rise to four distinct Dicer groups (DCL-1 to DCL-4) after 
the divergence of single-celled green algae but before or 
around the divergence of moss from higher plants (fig. 4). 
This result mirrors what we observed for animal Dicer — rapid 
diversification of the Dicer gene family coincident with 
the origins of multicellularity — suggesting that diversifica- 
tion of RNAi pathways may have played important and 
parallel roles in the development of multicellular animals 
and plants. 



As with our analysis of animal Dicers, our plant Dicer 
phylogeny was robust to different analysis methods, align- 
ment strategies, support measures, and evolutionary models 
(supplementary table S1 and fig. S13, Supplementary Material 
online). The same tree topology was recovered using 
heterogeneous evolutionary models (supplementary figs. 
S14 and S15, Supplementary Material online) and alignment 
of the DCL "catalytic core" (supplementary fig. S16, 
Supplementary Material online), suggesting that our plant 
phylogeny is reliable. 

Although DCL-1 and DCL-3 have clear moss orthologs in 
all our analyses, suggesting this gene duplication occurred 
before the divergence of moss from higher plants, the precise 
timing of the DCL-2/DCL-4 duplication was unresolved in 
some trees (fig. 4). To examine the DCL-2/4 duplication 
in more detail, we aligned additional Dicer homologs from 
a variety of plant genomes and reconstructed the plant DCL 
phylogeny using this expanded dataset. 

The resulting tree identified a moss DCL-4 ortholog and a 
Selaginella DCL-2 ortholog, both with fairly strong statistical 
support (supplementary fig. S17, Supplementary Material on- 
line), consistent with what we found in our phylogeny built 
from the DCL "catalytic core" (supplementary fig. S16, 
Supplementary Material online). Although we remain cau- 
tious when interpreting these findings — as the timing of 
the DCL-2/4 duplication was unresolved in some trees — 
these results suggest that the DCL-2/4 duplication may 
have also occurred before moss diverged from higher 
plants, in which case all four DCLs would have originated 
very early in plant evolution, coincident with the origin of 
multicellularity. 

We observed some changes in the Dicer domain architec- 
ture across the plant phylogeny (fig. 4) — similar to what 
we observed in animals — including lineage-specific losses of 
functional DEAD/Helicase domains in monocot and dicot 
DCL-1, which is known to function primarily in miRNA- 
based gene regulation (Welker et al. 2011). The major differ- 
ence between plant and animal Dicer domain architecture 
appears to be the presence of twin C-terminal dsrm domains 
in many plant Dicers; this architecture was not observed 
in the animal proteins. Twin dsrm domains were found 
throughout DCL-4 proteins and in DCL-1 s. 

The functional consequences of these changes in 
C-terminal dsrm architecture are unknown. Dsrm domains 
typically function in coordinating the hand-off of the RNA 
template from Dicer to an AGO protein, either through pro- 
tein-RNA or protein-protein interactions (Parker et al. 2008). 
This suggests that differences in C-terminal domain architec- 
ture may play a role in defining specific RNAi pathways by 
determining which downstream partners a specific Dicer 
interacts with (Marques et al. 2010). 

In contrast to the sparse protein-coding adaptation we 
observed for animal Dicers, we found strong evidence for 
recurrent adaptation targeting the PAZ domain of plant 
DCL-4 throughout monocots and dicots (fig. 4). There was 
little evidence for adaptation in other functional domains or 
other DCL groups, suggesting that adaptation has specifically 
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Fig. 4. Dicer duplicated early in plant evolution. We inferred the plant Dicer family phylogeny using maximum likelihood and Bayesian methods. 
Support is shown as SH-like aLRT scores / bootstrap proportions / Bayesian posterior probabilities (supplementary figs. S13-S17, Supplementary 
Material online). We additionally display the inferred domain architecture for each sequence. Colored circles on each branch indicate significant support 
(P < 0.05 after correcting for multiple tests) for adaptive protein-coding changes in each functional domain, inferred using branch-sites analysis (Zhang 
et al. 2005). Species names are colored by taxonomic group. 
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driven changes in antiviral DCL-4 PAZ function throughout 
plant evolution. 

The Dicer PAZ domain anchors the end of the RNA mol- 
ecule to facilitate cutting by twin RNase3 domains. As is the 
case with many similar RNA-binding domains, Dicer PAZ 
contains a large, positively charged pocket that binds the 
RNA-end structure primarily through electrostatic inter- 
actions (Ma et al. 2004). As expected, we found that this 
pocket was highly conserved across Dicers, including both 
plant and animal proteins (fig. 5). However, the plant 
DCL-4 PAZ domain showed extreme variation in electrostatic 
distribution across the RNA-binding pocket, even among very 
closely related plants. 

Particularly in monocots, there appears to be a complete 
reversal of the electrostatic charge of the DCL-4 PAZ RNA 
binding pocket to a primarily negative charge. This radical 
change in pocket charge is completely novel among Dicer 
PAZ domains and is uncharacteristic of PAZ domains in 
general (Lingel et al. 2004; Wang, Juranek, et al. 2009). 
Although the functional consequences of these changes in 
binding-pocket electrostatics are unknown, such extreme 
changes are highly likely to affect RNA binding, suggesting 
that DCL-4 may exhibit radically different RNA-binding prop- 
erties compared with other Dicer proteins. 

Given the importance of plant DCL-4 for antiviral immun- 
ity (Bouche et al. 2006; Deleris et al. 2006; Diaz-Pendon et al. 
2007) and the probable role of the PAZ domain in recognizing 
specific RNA motifs (Lingel et al. 2004; Ma et al. 2004; 
Wang, Juranek, et al. 2009), our results suggest that the 
DCL-4 PAZ domain is likely involved in a long-term evolu- 
tionary arms race with viral factors, leading to rapid vari- 
ation in RNA-binding properties. When we identified 
specific adaptive substitutions in various DCL-4 homologs, 
we found that many of them clustered around the 
RNA-binding pocket, with a particularly high concentration 
found in the RNA-binding loop between helices 1 and 
2 (fig. 6). One position in the middle of this loop is highly 
variable across DCL-4s and has undergone adaptation in 
multiple plant lineages; this position is extremely conserved 
across other Dicer proteins (supplementary fig. S18, 
Supplementary Material online) and forms a critical 
RNA-contact in Giardia Dicer (Simon et al. 2011). The radical 
difference in both conservation and amino acid physiochem- 
ical properties at this critical RNA-contact residue between 
DCL-4 and other Dicers further suggests that DCL-4 
RNA-binding properties are likely to be distinct from those 
of other Dicers, although the precise nature of these func- 
tional differences is unknown. 

Evolution of miRNA Recognition 
Recent analyses suggest that miRNA genes may have 
originated independently in animals and plants, suggesting 
that miRNA-based gene regulation may have arisen at least 
twice in eukaryote evolution (Murphy et al. 2008; Lee et al. 
2010; Axtell et al. 2011). Dicer recognizes the 3' -end of 
miRNAs via specific contact residues within the PAZ 
domain (Zhang et al. 2002, 2004). Efficient processing of 
miRNAs in humans requires additional contacts with the 



5' -end of the miRNA molecule, facilitated by conserved 
residues within the C-terminal part of the PAZ domain as 
well as an N-terminal PAZ extension. Recent analyses suggest 
that human miRNA processing relies primarily on 5' miRNA 
recognition (Park et al. 2011). 

We examined the patterns of sequence evolution and 
intron/exon structure around key miRNA-contact residues 
in animal and plant Dicers to better understand how 
miRNA processing evolved. We found that a number of 3'- 
contacting residues and surrounding PAZ-domain sequence 
were highly conserved across not only all animal Dicerls but 
also all plant DCL-1 proteins, suggesting that the biochemical 
basis for 3' miRNA recognition was present in the ancestral 
Dicer before the animal/plant split (fig. 7). Maximum- 
likelihood ancestral reconstruction confirmed this conclusion, 
as all four ^-contacting residues were inferred in animal, 
plant, and animal/plant/fungi ancestral sequences with 
strong support (supplementary table S4, Supplementary 
Material online). These contact residues and surrounding se- 
quence were not conserved in animal Dicer2s or plant DCLs 
2-4, suggesting that the mode of RNA recognition and/or 
binding differs between miRNA-processing and other Dicers. 
Interestingly, some of the Dicer homologs falling in the unre- 
solved base of the animal tree had this conserved Dicer1-like 
signature, suggesting that they may have retained the ability 
to bind the 3'-end of miRNAs (fig. 7). 

Sequence alignment of the 5' miRNA-binding pocket 
revealed that while S'-binding residues within the PAZ 
domain (R996 and R1003) are conserved across eukaryote 
Dicers (but not other PAZ-containing proteins, supplemen- 
tary fig. S19, Supplementary Material online), key RNA- 
contact residues within the N-terminal extension (R778, 
R780, and R811) are conserved only among bilaterian 
Dicerl proteins and absent from other animal and plant 
Dicers (fig. 7). Gene modeling identified a bilateria-specific 
exon encoding the entirety of the S'-binding N-terminal 
extension sequence that is not present in other animal or 
plant Dicers, suggesting that this 5'-binding pocket is a bila- 
terian innovation. 

Although all plant Dicer homologs lack the bilaterian 
N-terminal extension of the PAZ domain required for effi- 
cient 5' miRNA binding, we identified a similar but 
non-homologous N-terminal extension of the PAZ domain 
conserved throughout plant DCL-1 but absent from other 
DCLs (fig. 7). This DCL-1 insertion does not align to the 
animal 5' miRNA-binding pocket and has different intron/ 
exon boundaries, supporting the conclusion that it is not 
homologous to the animal pocket but represents a unique 
plant DCL-1 sequence. 

Both the animal 5'-pocket and the plant DCL-1 insertion 
contain a large number of positively charged residues 
(Arg and Lys), which in animals have been shown to be im- 
portant for miRNA-binding (Park et al. 2011). Although func- 
tional analyses will be required to confirm the role of 
the DCL-1 insertion, evolutionary analysis suggests that 
the plant DCL-1 insertion may be a novel motif supporting 
efficient DCL-1 -based gene regulation by facilitating 5' 
miRNA-binding. 
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dicots 



monocots 



Fig. 5. Adaptive substitutions altered the electrostatic distribution across the RNA-binding pocket of plant DCL-4 PAZ. We inferred the 3D structure 
and electrostatic distribution of plant and animal Dicer PAZ domains (see Materials and Methods). We plot the electrostatic distribution (kT/e) across 
the RNA-binding pocket (yellow-dotted outline) of each protein. 



635 



Mukherjee et al. • doi:10.1093/molbev/mss263 



MBE 




Adaptive substitutions 

• O. sativa - S. bicotor 

• O. sativa - Z. mays 

• A. thaliana - B. rapa 

RNA-contact residues 



Qryza sativa 
Zea mays 
Sorghum btcotor 
Brassica rapa 
Arabidopsis thatiana 



Qryza sativa 
Zea mays 
Sorghum bicolor 
Brassha rapa 
Arabidopsis thaliana 



m m 

ILLDGI F 5 K T 
L LGGTYNKA 
L LGGAYKKA 
L ANGCW5 I D 
L A NGCWN I D 



* * ★# • ** • * 

VVGSVV F SP HNN I F F F vIg I lIe I l AW H| s l A l Y A 

V IGSLVYTPHTDL FF F V ©|A I L D G T N AK S1LNGA Iy A 

V IGSLVYTPHTDLFF F v| S I l_HGT N AK SflL NG A|Y E 

VKNSVVFATHKKQFYFVSDICHGRNGFSA 
VKN5 LV FTTYSKQ FY FVAD 



I 



H F K | 



HFKER 



I C HG R NG F S A I RKS ITETHL 

ichgrngfspvkesItkshv 



E S 



★* 


★ 


* 




F R - - - 


I E L S H P 


EQP 




I H 


I K L SHP 


EQP 


ill 


1h- - - 


1 EL SHP 


EQP 




In s y| 


VKLKHP 


LQP 




|k ly| 


v|l k hp 


AQP LLR 




LEH 


FV 


M EH 


F V 


M E H 


F V 


FEY 


F 1 


D E Y 


F 1 







E i P P 
1 1 P P 

Y P 





1 G F 




1 G F 


L K 1 


T G F 




K G L 




KGL 



IKD. 

;kdi 

t K Dl 

; k d 



MGSSLSLLPS L M Y | 

MGSSLSLLPS LMyI 




GSSLSLLP S VMH 
GSSLSLLPS 



Fig. 6. Adaptive substitutions in plant DCL-4 PAZ domain are likely to affect RNA binding. We inferred adaptive protein-coding substitutions in 
monocot and dicot DCL-4 PAZ domains using branch-sites models (see Materials and Methods). We plot adaptive substitutions (posterior probabil- 
ity > 0.95) along the protein structure (top) and multiple sequence alignment (bottom). RNA-contacting residues identified in Giardia Dicer PAZ 
(Simon et al. 2011) and AGO PAZ (Wang, Juranek, et al. 2009) are indicated below the sequence alignment. 



Discussion 

Our results suggest a model of Dicer evolution in which an 
ancestral eukaryote Dicer duplicated independently and early 
in animal and plant lineages, coincident with the origins 
of multicellularity. These early Dicer family expansions may 
have coincided with duplications of other proteins in the 
RNAi pathway (such as AGO, supplementary fig. S12, 
Supplementary Material online) as well as with the independ- 
ent origins of miRNA genes in animals and plants (Murphy 
et al. 2008; Lee et al. 2010; Axtell et al. 2011), suggesting that 
the requirement for complex gene regulation to support 
multicellular organisms may have independently recruited 
RNAi pathways and driven their elaboration in animals and 
plants. 

Our results suggest that the major functional repertoire of 
the Dicer family was likely present in early multicellular ani- 
mals and plants, with subsequent functional refinement and 
gene loss events "fine-tuning" RNAi function in various lin- 
eages. Specifically, bilaterian Dicerl appears to have evolved 
a novel 5' miRNA binding mechanism, which increased 
miRNA-processing efficiency, a trait that may have arisen in- 
dependently in plant DCL-1. We also observed parallel losses 
of the DEAD/Helicase domain of Dicerl and DCL-1 in various 
protostome and plant lineages, which likely focused these 
proteins' functions on miRNA processing. 



Although the importance of RNAi for antiviral immunity 
in vertebrates is unclear (Cullen 2006), Dicer2 and DCL-4 are 
known to play critical antiviral roles in invertebrate animals 
and plants, respectively (Hamilton and Baulcombe 1999; 
Zamore et al. 2000; Aliyari and Ding 2009). Our model of 
animal Dicer evolution suggests that Dicer2's antiviral role 
probably arose very early in animal evolution, with Dicer2 
being subsequently lost from deuterostomes and nematodes. 
This raises the possibility that alternative antiviral strategies 
may have compensated for loss of RNAi-based antiviral im- 
munity in these lineages. 

Although we detected no evidence for long-branch attrac- 
tion or other phylogenetic artifacts from any of the major 
factors known to cause such errors, there is always the pos- 
sibility that the inferred phylogeny is not the correct tree. In 
the case of the animal Dicer phylogeny, the presence of mul- 
tiple long Dicer2 branches falling at the base of the tree is 
consistent with a long-branch attraction hypothesis in which 
fast-evolving Dicer2 duplicates — perhaps arising in multiple 
lineages — artifactually group at the base of the animal tree. 
More extensive taxon sampling of basal animal lineages will 
help resolve this issue. 

Perhaps, the most striking result from our study is the 
extremely strong signature of recurrent adaptation through- 
out the evolution of the plant DCL-4 PAZ domain and the 
radical variability in its RNA-binding pocket compared with 
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Fig. 7. Examination of Dicer sequence alignment reveals bilaterian- and DCL-1 -specific N-terminal PAZ extensions likely to contribute to 5 r 
RNA-binding. We display parts of aligned Dicer sequences from representative taxa. Protein domain architecture is show in the middle, with the 
PAZ domain indicated in yellow. Specific residues previously shown to contribute to 3' and 5 r RNA-binding are indicated at the top of each alignment 
(Lingel et al. 2004; Wang, Juranek, et al. 2009; Park et al. 2011; Simon et al. 2011). Alignment of conserved N-terminal PAZ extensions from bilateria and 
plant DCL-1 are indicated in pink and orange, respectively. 
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those of other Dicers. The PAZ domain anchors the end of 
the RNA molecule and may play an important role in recog- 
nizing unwanted RNAs (Simon et al. 2011). Most of the adap- 
tive protein-coding substitutions we identified in DCL-4 PAZ 
are located around the RNA-binding pocket, and some affect 
key RNA-contact residues, suggesting that changes in viral 
RNA structure may have driven adaptation in DCL-4 PAZ. 
Alternatively, several viral suppressors of RNAi (VSRs, viral 
proteins functioning to disrupt the host's RNAi-based 
immune response) are known to antagonize the PAZ 
domain and may have driven its rapid evolution in plants 
(Baumberger et al. 2007; Singh et al. 2009; Hamera et al. 2012). 
Uncovering the precise interactions driving this clear evolu- 
tionary arms race will help determine the general biochemical 
forces governing host-pathogen interactions. 

RNAi biochemistry operates on the elegantly simple 
principle of using an RNA template in coordination with 
enzymatic proteins to target and regulate complementary 
RNA molecules. Organisms have exploited this simple 
principle to support multiple functional roles, such as gene 
regulation, antiviral immunity and preserving genome integ- 
rity. Our analyses have revealed the RNAi system to be highly 
dynamic across evolutionary timescales, with parallel gains 
and losses of entire functional pathways, long-term evolution- 
ary arms races, and changes in protein domain architecture 
and sequence driving functional divergence. Given this 
dynamic backdrop, we expect RNAi to function differently 
in different organisms; future comparative-functional studies 
are expected to further illuminate the details of how evolu- 
tion has shaped RNAi function. 

Materials and Methods 

Gene Identification 

We used confirmed functional Dicer genes as queries to 
identify new Dicer homologs from fully sequenced genomes 
using TBLASTN, with a liberal e-value cutoff of 1.0 to identify 
all potential homologs (Altschul et al. 1997). Sequences 
recovered from one round of TBLASTN were recursively 
used as queries until no further sequences were detected. 
We additionally performed TBLASTN searches using individ- 
ual functional domains as query sequences, to identify poten- 
tial Dicer homologs with different domain architectures. 

We used homology-based gene prediction based on 
hidden Markov models (HMMs) to identify the complete 
open reading frame of each potential Dicer homolog. 
We excised ±2,000 nucleotides of genomic sequence sur- 
rounding each TBLASTN hit and calculated log-odds scores 
for coding exons using FGENESH+ with a log-odds cutoff 
of 50.0 to identify strongly supported exons (Salamov and 
Solovyev 2000). Potential exon mis-predictions were manually 
refined using dot-plot analysis and 3-frame translations 
of the genomic sequence. The closest known homolog 
recovered from PSI-BLAST search of Gen bank was used as a 
guide for homology-based gene prediction. Coding exons 
were re-predicted, if necessary, following multiple sequence 
alignment of the resulting protein sequences (see below), and 
non-homologous sequences were removed. 



Protein Domain Identification 

Protein functional domains were identified by sequence 
search of the PFam database (Punta et al. 2011), using an 
e-value cutoff of 1.0e -5 to identify strongly supported 
domains. PFam hits were verified via sequence searches of 
the SMART (Letunic et al. 2012) and NCBI conserved domain 
databases (March ler-Bauer et al. 2011). We also manually 
verified loss of functional DEAD domains by examining 
the highly conserved "DECH" motif in each sequence. 
We identified Dicer 5' RNA-binding pockets by aligning 
Dicer PAZ domains and relevant N-terminal sequences and 
manually examining these alignments for the presence or 
absence of experimentally determined RNA-binding residues 
(Park et al. 2011). 

Sequence Alignment 

Sequences were aligned using Muscle V3.8.31 (Edgar 2004) 
(default parameters), Mafft v6.850b (Katoh et al. 2002) 
(1,000 iterations and -genafpair option) and Kalign v2.03 
(Lassmann and Sonnhammer 2005) (default parameters). 
All alignments were manually refined to remove obvious 
alignment errors and preserve alignment of identified func- 
tional domains. Resulting alignments were processed using 
Gblocks v0.91b (Castresana 2000) to remove ambiguously 
aligned regions prior to phylogenetic analysis, allowing gaps 
in 1 /2 of the sequences. Additionally, we produced an elision 
alignment by concatenating all three individual Gblocks- 
processed alignments (Wheeler et al. 1995). In addition to 
full-length and elision alignments, we produced individual- 
domain alignments consisting of only the DEAD/Helicase 
domains and only the Dicer "catalytic core" PAZ+RNase3 
domains. These individual-domain alignments were not 
processed using Gblocks. 

Phylogenetic Analysis 

Maximum likelihood trees were inferred using PhyML v3.0 
(Guindon and Gascuel 2003; Guindon et al. 2010), with the 
best-fit evolutionary model being selected using ProtTest v3.0 
(Darriba et al. 2011). Clade support was calculated using 
SH-like approximate likelihood ratio tests (Anisimova et al. 
2011) and bootstrap proportions (500 replicates). 

Bayesian phylogenies were reconstructed using MrBayes 
v3.1.2 (Ronquist et al. 2012), assuming the default priors. 
We integrated over all available substitution models and 
assumed a discrete-gamma model of among-site rate vari- 
ation. We evaluated chain convergence by performing two 
independent runs and terminating the analysis when the 
average standard deviation in clade probabilities between 
the two runs fell below 0.01. 

We also reconstructed phylogenetic trees using two 
approaches that have been shown to improve phylogenetic 
accuracy. First, we used a Bayesian approach implemented in 
PhyloBayes v2.3 that incorporates variation in amino acid 
frequencies across the sequence alignment (Lartillot et al. 
2009). Second, we used a maximum-likelihood method that 
incorporates across-site heterogeneity in branch lengths 
(Kolaczkowski and Thornton 2008). For this approach, the 
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best-fit number of branch-length categories was identified 
using the Akaike information criterion (AIC). 

Protein Adaptation 

We evaluated protein-coding adaptation across the phyl- 
ogeny using a branch-sites test implemented in PAML v4.5 
(Zhang et al. 2005; Yang 2007). We used a four-parameter 
branch-sites model: a proportion of sites in the alignment (p 0 ) 
are assumed to be under purifying selection, evolving at non- 
synonymous/synonymous rate 0 < w 0 < 1. A second class of 

sites are neutral (n/ q = 1), and a third class of sites 
(p 2 =[1 — Po — Pi]po/[Po + Pi]) evolve under positive selec- 
tion (yv 2 > 1). The null hypothesis constrains w 2 - 1. 

For each branch on the phylogeny, we calculated the like- 
lihood ratio of the positive-selection model (w 2 > 1) versus 
the constrained neutral model (u/ 2 = 1) and estimated 
P values using the chi-square distribution with one degree 
of freedom (Zhang et al. 2005). We corrected for multiple 
testing using a Bonferroni correction. We additionally con- 
firmed any positive results by re-testing after removing 
fast-evolving taxa, which might cause false positives. We 
report only those results that were positively identified 
(P < 0.05) in both analyses. 

Individual codons were classified as adaptive or 
non-adaptive using Bayes empirical Bayes (BEB) posterior 
probabilities calculated by PAML. For each codon, we com- 
pare the probability that the codon evolved with w 2 > 1 in a 
given lineage to the sum of the probabilities that the codon 
evolved with 0 < w 0 < 1, ma, = 1 or w 2 > 1, each weighted by 
the estimated proportion of sites evolving under that cat- 
egory. Parameter uncertainty is incorporated by integrating 
over diffuse prior distributions (Zhang et al. 2005). Sites with 
BEB posterior probability >0.95 were considered adaptive 
protein-coding changes. 

Structural Modeling 

We identified homologous X-ray structures via sequence 
search of the Protein Data Bank (Bernstein et al. 1978). 
Sequences were aligned to the structural template manually 
following initial multiple-sequence alignment by Muscle. We 
built five structural models of each sequence using Modeller 
v9.10 (Sali and Blundell 1993) and report results obtained 
using the best model, selected based on DOPE score (Shen 
and Sali 2006). 

Structural models were processed using PROPKA and 
PDB2PQR to determine residue side-chain pKas, optimize 
the structure for favorable hydrogen bonding and calculate 
charge and radius parameters from electrostatic force fields 
(Dolinsky et al. 2004; Rostkowski et al. 2011). Electrostatic 
surface potentials were estimated from processed structures 
using APBS (Baker et al. 2001) and projected onto the mo- 
lecular surface for visualization. 

Supplementary Material 

Supplementary Excel spreadsheet, tables S1-S4, and figures 
S1-S19 are available at Molecular Biology and Evolution online 
(http://www.mbe.oxfordjournals.org/). 
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